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THE BLACK HOLE MASS, STELLAR M/L, AND DARK HALO IN M87 
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ABSTRACT 

We model the dynamical structure of M87 (NGC4486) using high spatial resolution long-slit observations of 
stellar light in the central regions, two-dimensional stellar light kinematics out to half of the effective radius, and 
globular cluster velocities out to 8 effective radii. We simultaneously fit for four parameters, black hole mass, 
dark halo core radius, dark halo circular velocity, and stellar mass-to-light ratio. We find a black hole mass of 
6.4(±0.5) x 10 9 Mq (the uncertainty is 68% confidence marginalized over the other parameters). The stellar 
M/Ly = 6.3 ± 0.8. The best-fitted dark halo core radius is 14 ±2 kpc, assuming a cored logarithmic potential. The 
best-fitted dark halo circular velocity is 715 ± 15 km s" 1 . Our black hole mass is over a factor of 2 larger than 
previous stellar dynamical measures, and our derived stellar M/L ratio is 2 times lower than previous dynamical 
measures. When we do not include a dark halo, we measure a black hole mass and stellar M/L ratio that is 
consistent with previous measures, implying that the major difference is in the model assumptions. The stellar 
M/L ratio from our models is very similar to that derived from stellar population models of M87. The reason for 
the difference in the black hole mass is because we allow the M/L ratio to change with radius. The dark halo is 
degenerate with the stellar M/L ratio, which is subsequently degenerate with the black hole mass. We argue that 
dynamical models of galaxies that do not include the contribution from a dark halo may produce a biased result 
for the black hole mass. This bias is especially large for a galaxy with a shallow light profile such as M87, and 
may not be as severe in galaxies with steeper light profiles unless they have a large stellar population change with 
radius. 

Subject headings: galaxies: elliptical and lenticular, cD; galaxies: individual (M87, NGC4486); galaxies: 
kinematics and dynamics 



1. INTRODUCTION 

M87 has been the poster child for black hole studies due to 
its central jet, large central velocity dispersion, and having the 
largest measured black hole mass. Sargent et al. (1978) and 
Young et al. (1978) estimate a black hole mass of 5 x 10 9 M Q 
using the surface brightness and velocity dispersion profiles. 
Subsequently, a large number of studies based on the stellar dy- 
namics provide estimates from 1-3 x 10 9 Mq (see review of 
Kormendy & Richstone 1995). The latest estimate from Harms 
et al. (1994) use Hubble Space Telescope (HST) to provide a 
gas kinematic map and black hole estimate based on gas dy- 
namics. Macchetto et al. (1997) analyse additional HST gas 
kinematics to derive a black hole mass of 3.2(±0.9) x 10 9 M Q , 
which is considered the standard estimate for M87. 

As one of the nearest galaxies and one of the largest galax- 
ies, it is clear why M87 has had a long and important history 
of having it's black hole mass measured. In fact, M87 is often 
used as the anchor for the upper end of the black hole mass dis- 
tribution. There are significant consequences for understanding 
the value of the maximum black hole mass. It strongly effects 
the parameterization of the black hole mass correlations (e.g., 
Gultekin et al. 2009), the space density of the most massive 
black holes (e.g., Lauer et al. 2007), the comparison with mass 
derived from quasars which tend to be closer to 10 10 rather than 
10 9 Mq (e.g., Shields et al. 2003), among many others. 

A significant issue for stellar dynamical models of black hole 
in M87 is that the dark halo contribution has never been in- 
cluded in the modeling. The latest model of the gravitational 
potential come from Wu & Tremaine (2006), McLaughlin 
(1999), Kronawitter et al. (2000), Romanowsky & Kochanek 



(2001), and Churazov et al. (2008). Wu & Tremaine use the 
globular cluster kinematics to measure the enclosed mass as- 
suming a two-integral distribution function. McLaughlin use a 
isotropic model for M87 based on the globular cluster kinemat- 
ics and do not include a black hole. Both Kronawitter et al. and 
Romanowsky & Kochanek (2001) assume a spherical distribu- 
tion. Churazov et al. use the gravitational potential as derived 
from X-ray gas and compare to that derived from stellar dynam- 
ics. These studies find moderate agreement for the gravitational 
potential at large radii from the two techniques. Our goal in this 
paper is to understand possible degeneracies between measur- 
ing the central black hole and the dark halo simultaneously. We 
use orbit-based models which do not limit the form of the al- 
lowed velocity anisotropies assuming axisymmetry. We focus 
on M87 since the quality of the observational data for M87 has 
dramatically increased, and it is worthwhile to analyze it again 
with the most up-to-date data and dynamical models. 

Our main result is that we find a black hole mass that is over 
2 times larger when including a dark halo compared to models 
that do not contain a dark halo. This result is both generic in 
that it could be a potential concern for many of the black hole 
studies to date, and specific in that it applies the most up-to-date 
modeling to M87. Humphrey et al. (2008) model NGC 4649, 
a galaxy similar in mass to M87, using X-ray gas emission and 
find a black hole that is 2 times larger than the previous stellar 
dynamical measure of Gebhardt et al. (2003). Thus, there is a 
concern that black hole masses at the upper end of the distribu- 
tion may be biased. 

We assume a distance to M87 of 17.9 Mpc. 
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2. DATA 

In the central regions of M87, we rely on the stellar light for 
both photometry and kinematics. For the outer region, the stel- 
lar kinematics only extend to about 40". Therefore, at large 
radii we rely on radial velocities of individual globular clusters 
for the kinematics, also including the number density profile 
for the clusters. Unfortunately, the stellar data ends around 40" 
and the high signal-to-noise globular cluster data starts around 
150". Thus, there is a radial region where we presently have no 
kinematic data. 
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FIG. 1. — The volume density radial profile for stars. The black and red 
solid line represent the luminosity density along the major and minor axis, 
respectively. These come from deprojection of the surface brightness profile 
presented in Kormendy et al. (2009). The dashed line is the deprojection of 
the surface number density for globular clusters as measured by McLaughlin 
(1999). The globular cluster profile is normalized, arbitrarily, to match the 
stellar profile in order to compare the shapes. The actual volume density of the 
clusters is much lower and their contribution to the stellar mass is irrelevant. 
Furthermore, the cluster profile within 20" is extrapolated since there is little 
information to measure a surface brightness profile. Both profiles are used in 
the dynamical models for each of the components. 

2.1. Stellar Density 

The stellar surface density profile comes from Kormendy et 
al. (2009), which is a combination of HST data from Lauer 
et al. (1992) and various ground-based observations. It ex- 
tends from 0.02" to 2400" (5 orders of magnitude). The sur- 
face brightness profile in the V-band, ellipticity, and position 
angle versus radius are shown in Kormendy et al. (2009). The 
ellipticity is near zero for most of the radial range, but then in- 
creases at both small and large radii. The increase at large radii 
has been known for some time (e.g., Weil et al. 1997), and is 
associated with the large rotation seen there (Kissler-Patig & 
Gebhardt 1998). The increase at small radii has been noted by 
Lauer et al. (1992) and appears to be real (as opposed to noise). 
In fact, the axisymmetric deprojection of the surface brightness 
to the luminosity density requires the increase at small radii. 
We deprojected the surface brightness using Magorrian (1999). 
This deprojection allows for a radial variation of the observed 
ellipticity. Figure 1 plots the luminosity density along the ma- 



jor and minor axis versus radius. The dynamical models use 
this non-parametric model of the luminosity density. From this 
figure one can see the increase in the ellipticity in the central 
region and in the outer region (the minor axis is lower than the 
major axis). From about 2-70" M87 is nearly round in projec- 
tion. The full radial range used in the dynamical models spans 
10 7 in the luminosity density. 

2.2. Stellar Kinematics 

There have been many stellar kinematic studies of M87 and 
we rely on the two best studies for our models. These are the 
kinematics from SAURON (Emsellem et al. 2004) and from 
van der Marel (1994). 

The SAURON Team provide two-dimensional kinematic 
coverage of M87 out to about 40". The kinematics are parame- 
terized with the first four coefficients of a Gauss-Hermite poly- 
nomial expansion, with uncertainties. The two-dimensional 
coverage of SAURON provides the best possible coverage for 
determining the stellar orbital structure and mass-to-light pro- 
file. Our dynamical models fit the line-of-sight velocity dis- 
tribution (LOSVD) directly as opposed to fitting to moments 
or Gauss-Hermite coefficients. Thus, we have to convert the 
SAURON velocity moments into LOSVDs; we do this by cre- 
ating Monte Carlo simulations based on the uncertainties of 
the Gauss-Hermite coefficients. This procedure assumes that 
the moments are uncorrected (further discussed in Section 4). 
The spatial resolution of the dynamical models is coarser than 
the SAURON data in the outer regions and similar in the in- 
ner regions (cf. Sec. I3.21 i. Where necessary, we rebin the 
SAURON data by averaging weighted by the uncertainties over 
all SAURON LOSVDs in one model bin. We do not average the 
Gauss-Hermite coefficients. Ideally, we would not apply any 
spatial binning and use each ground-based element indepen- 
dently (even though the kinematic data was already binned orig- 
inally), but this fine binning would take a significant amount of 
computing resources as to make the study unfeasible. We dis- 
cuss the binning further in Sec 13. 21 

The kinematic data with the best spatial resolution come 
from van der Marel (1994). They have 0.6" FWHM seeing 
with a 1" slit. The data extend to 10". Van der Marel pro- 
vides the first four Gauss-Hermite moments and we transform 
to LOSVDS as we did with the SAURON data. 

Another important dataset comes from Sembach & Tonry 
(1996), who obtain kinematics in M87 out to a radius of 130". 
We do not use their kinematics mainly since they only pro- 
vide velocity and velocity dispersion, without higher order mo- 
ments. Since the SAURON data are such higher S/N and pro- 
vide the first 4 moments of the LOSVD and since our modeling 
code fits LOSVDs (and not moments), in order to include the 
Sembach & Tonry data we would have to mock up the higher 
order moments, which would lead to additional uncertainty. 
Furthermore, there appears to be an offset in the velocity dis- 
persion from Sembach & Tonry compared to other kinematic 
studies as discussed in Romanowsky & Kochanek (2001). Fu- 
ture datasets at these radii in M87 using the full LOSVD will 
be very important for additional analysis. 

There are many other kinematic datasets, but the S/N and 
areal coverage of the SAURON data is superior to all other stud- 
ies that it will make little difference, if any, to the dynamical 
models. Given there are often systematic differences between 
various datasets, it is better to use just one dataset. We include 
the van der Marel data since these add significantly to the re- 
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covery of the black hole mass. The two datasets are consistent 
with and show no biases between each other at the radii where 
they overlap. 




FIG. 2. — The grid lines represent the spatial binning used for the dynamical 
models. The points are the locations of the 278 globular cluster radial veloci- 
ties. The black and red points represent clusters on opposite sides of the minor 
axis, and the red points have had their velocity flipped relative to the systemic. 
The colored spatial bins are the binning used for those particular cluster. Clus- 
ters outside those highlighted bins are not used in the models. There are 275 
cluster velocities used for the analysis. 

2.3. Globular Cluster Data 

The surface density profile for the globular clusters come 
from McLaughlin (1999), who compile the profile based on var- 
ious datasets. We deproject this profile using a non-parametric 
spherical inversion as described in Gebhardt et al. (1996). In 
Figure 1, we include deprojected, normalized number density 
profile of the globular clusters; we use an arbitrary scaling for 
the cluster profile in order to compare the profile shapes of the 
clusters and stars. The total mass of the clusters is irrelevant 
compared to the stars. The dynamical models only use the slope 
of the luminosity density in order to constrain the enclosed 
mass. At large radii the globular cluster profile is shallower 
than the stellar profile, thus it is important to use the cluster 
profile in the dynamical models. Within 20" there is very little 
information on the globular cluster profile and we must use an 
extrapolation there. The dashed line in Figure 1 shows what 
we used for the models. We have also tried 2 other extrapola- 
tions: one assumes a flat profile inside of 10" and the other has 
a slope inside of 20" that is equal to the stellar density. These 
two profiles give nearly identical results as the adopted profile. 

Globular cluster velocities are reported in Cohen (2000), Co- 
hen & Ryzhov (1997) and Hanes et al. (2001). These are com- 
piled in Cote et al. (2001). We use the same cut as defined in 
Cote et al. to remove foreground and background contamina- 
tion, based on velocity, color, and magnitude. These cuts result 
in 278 velocities. 

We include the individual velocities in the dynamical models 
by measuring the LOSVD in different spatial regions. We di- 
vide the 278 velocities into 1 1 spatial bins. These spatial bins 



are shown in Figure 2, where we plot the model grids. Due 
to the nature of the spatial binning for the dynamical models, 
some bins will have too few clusters to be useful. At radii be- 
tween 40 and 140" the density of clusters is low enough that we 
include larger bins in angle and radius than for the outer bins. 
Of the 278 velocities, we then use 275. 

Due to the uncertain number density profile for the globular 
clusters in the central regions, we have run models including 
and excluding clusters inside of 170". There are 68 cluster ve- 
locities in this range. Inside of 20", there are so few clusters 
known (since the galaxy light begins to dominate over the clus- 
ter light), that we already have to extrapolate the number den- 
sity profile. Inside of 100", the number density is better known, 
but since the light profile for M87 is so shallow, deprojection er- 
rors can lead to a large effect. There are no significant changes 
to the besst-fitted values between including or excluding the 
central cluster velocities, but the uncertainties improve when 
we use more data. The largest difference, which we discuss 
later, is that the dark halo scale radius is much better measured. 

From the approximately 25 velocities per spatial bin, we 
measure an LOSVD and uncertainties. The LOSVDs use 
an adaptive kernel density estimate adapted from Silverman 
(1986), and explained in detail in Gebhardt et al. (1996). The 
68% confidence bands for the LOSVD come from bootstrap re- 
samplings of the data. Another option would be to calculate 
moments from the individual velocities and compare the mod- 
els to the moments. Since the dynamical models fit LOSVDs 
directly, we choose to use the measured LOSVDs. The mo- 
ments derived from the measured LOSVDs compare very well 
with the moments derived from the individual velocities. 

Since we are using an axisymmetric model, the velocities 
are folded along the major axis and flipped about the minor 
axis (the velocity relative to the systemic velocity switches sign 
when flipped to the opposite side of the minor axis), such that 
all data are placed in one quadrant. This folding preserves any 
rotation or angular differences at a given radius. Rotation for 
the globular clusters is significant at large radii (Kissler-Patig 
& Gebhardt 1998, Cohen 2000, Cote et al. 2001) and needs to 
be included for a proper dynamical model. 

3. DYNAMICAL MODELS 

The orbit-based models used here are described in Geb- 
hardt et al. (2000, 2003), Siopis et al. (2009) and Thomas 
et al. (2004, 2005), and are based on the idea presented in 
Schwarzschild (1979). Similar orbit-based models are pre- 
sented, for example, in Rix et al. (1997), van der Marel et al. 
(1998), Cretton et al. 1999, Valluri et al. (2004). The modeling 
procedure consists of five steps which are further detailed be- 
low: (i) the surface-brightness distribution is deprojected into 
a 3d luminosity profile; (ii) a trial gravitational potential is cal- 
culated including - in our case here - the contribution from 
the stars, a dark matter halo and a black hole; (iii) a large set 
of orbits is run in the trial potential; (iv) an orbit superposi- 
tion model of the galaxy in the trial potential is constructed that 
matches the kinematical data as good as possible subject to the 
constraints given by the 3d light profile and the trial potential; 
(v) steps (i)-(iv) are repeated for a large variety of trial mass 
profiles. A x 2 -analysis then determines the best-fit model and 
its uncertainties. 

These models have undergone extensive tests for both recov- 
ery of the dark halo and the black hole mass. Thomas et al. 
(2005) describe the ability to recover the dark halo properties. 
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Gebhardt (2004) and Siopis et al. (2009) describe the accuracy 
on the recovery of the black hole mass and the stellar M/L ratio. 
The orbit-based models show no significant bias in the recovery 
of the parameters for axisymmetric systems. 

3.1. Model assumptions about the mass distribution in M87 

We calculate two different sets of dynamical models for M87 
based on different assumptions about its mass structure. 

In a first set of models we do not consider a dark matter halo 
explicitly. In fact, we assume that the mass-to-light ratio of the 
galaxy is spatially constant (except at the origin): 

p=?v+M BH 6(r), (1) 
with Mbh being the mass of the central supermassive black hole 
and v being the 3d light profile. If it would be true that there 
is no dark matter around M87 and if the mass-to-light ratio of 
the stars in M87 would be constant throughout the galaxy, then 
T would be our trial value for the stellar mass-to-light ratio. 
However, in fact we do expect dark matter around M87 and T 
from equation (fTJ is supposed to be generally larger than the 
true stellar mass-to-light ratio, because it will also account for 
the potential dark mass. 

In our second set of models we assume that the mass distri- 
bution p of M87 can be written as the sum of (1) the stellar mass 
density, (2) a dark matter halo and (3) a central black hole: 

p=Tv + p DM +M BH 5(r). (2) 



The dark halo is included as described in Thomas et al. (2005). 
We use both a logarithmic halo with a density profile given as: 
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and a NFW potential as given by: 
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(Navarro, Frenk & White 1996). The logarithmic halo features 
a flat central density core of size r c and an asymptotically con- 
stant circular velocity v c . The NFW profile diverges like r~ l 
towards the center and drops off with r oc r" 3 in the outer parts, 
steeper than the logarithmic halo. Concentration c, scale-radius 
r s and the virial radius r v of NFW halos are related via c = r v /r s . 

3.2. Modelling the stellar kinematical data 

For the dynamical models of the stellar kinematics we use 
a spatial binning with N r = 28 radial and N$ = 5 angular bins, 
respectively (cf. Gebhardt et al. 2000 for more details on the 
spatial grids). Figure 2 plots the radial and azimuthal plot of 
the spatial binning. The model bins are generally larger than 
the SAURON bins and, in that case, we average the SAURON 
LOSVDs as described in Sec. 12.21 After the spatial binning we 
have TYf rs =6+46 LOSVDs (van der Marel & SAURON). Each 
of these stellar LOSVDs £ stars is sampled by 7V vel = 19 velocity 
bins. 
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FIG. 3. — Velocity v, velocity dispersion a and Gauss-Hermite coefficients h$, as a function of radius. Points with errors bars: observational data from van 
der Marel (red; left column) and SAURON (grey; columns on the right-hand side; each column corresponds to one angular bin of the models; the mean lattitude is 
quoted on top); solid lines: best-fit model with dark matter halo; dotted lines: best-fit model without dark matter. Bottom row: \ 2 difference between the best-fit 
with dark halo (Xqm) an d without (Xsc)- Positive values in the bottom row indicate that the model with dark matter fits better. 
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For each trial potential (cf. equations Q] and [2]i we determine 
the weighted orbit superposition that has the lowest 



-^stars 




(5) 



Here, £™ od M is the model LOSVD in spatial bin j. During the 
X 2 -minimization, the orbit model is forced to reproduce the 3d 
luminosity profile v with machine precision. This ensures that 
the final model is self-consistent, i.e. that the generated dynam- 
ical model supports those orbits out of which it is constructed. 

When two orbits differ in their spatial shape only on a sub- 
bin scale, then they are equivalent with respect to the density 
constraints, but their projected kinematics may be different. 
Then, a slight change in the potential parameters can be com- 
pensated for by shifting light from one orbit to the other, with- 
out affecting the density constraints. That our model bins are 
larger than the SAURON data therefore likely increases the un- 
certainties of the derived mass parameters with respect to the 
optimum that could have been achieved given the observational 
data. However, the overall mass budget of an axisymmetric 
dynamical model is determined (1) by the density constraints 
and (2) by the total amount of kinetic energy along the line-of- 
sight (Thomas et al. 2007a). Since Fig. 3 shows that there is 
no significant variation of the velocity dispersion over any of 
our model bins, the total kinetic energy along the line-of-sight 
is conserved after the re-binning. Furthermore, Siopis et al. 
(2009) and Thomas et al. (2005) perform tests of recovery of 
black hole mass and dark halo profile for analytic models. Us- 
ing the same modelling code as used here and a similar binning 
scheme, they find accurate recovery of the input parameters. 
Thus, we do not expect that the chosen spatial binning used 
here introduces a significant bias to the mass recovery. 

3.3. Modelling the globular cluster kinematics 

The globular cluster system can have a different orbital struc- 
ture than the system of the stars, such that we need to model the 
globular cluster velocities separately from the stars. However, 
the models for the GC system are very similar to the models 
of the stellar kinematics. In particular, we use the same mass 
models (cf. equation [2] we do not attempt to model the GC 
kinematics without a dark matter halo). The only difference is 
that the weighted orbit superposition in each trial potential is 
determined by minimizing 
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whe re C GC are the N GC = 11 globular cluster LOS VDs from 
Sec. 12.31 and «gc is the deprojected number-density profile of 
the GCs. When minimizing the above Xgc we f° rce the mod- 
els to reproduce «gc with machine precision. 

In contrast to the stars, we treat the GCs of M87 as a sys- 
tem of massless test particles. Therefore, in principle we do 
not need to reproduce «gc with machine precision, as any mis- 
match in «gc wouldn't have any effect on the potential. How- 
ever, we have tried three profiles for «gc with different central 
slopes and find nearly identical results. As discussed earlier, 
the more important aspect is whether we include the globular 
cluster velocities inside of 170" where the density of observed 
velocities begins to decrease. The analysis presented here use 
all of the central globular cluster velocites. However, models 




that exclude the central velocities show no change in the best- 
fitted values, but increase the uncertainties. The largest change 
is for the dark matter scale radius, since when including the cen- 
tral velocities, we find an increase in the uncertainty by about 
50%. 

3.4. Constraining the mass density in M87 

In total we fit for four sets of models: 1) we fit only for 
the black hole mass and stellar T including the SAURON and 
van der Marel data only (no globular cluster data and assum- 
ing no dark halo), which mimics what has been done for most 
black hole mass studies, 2) we fit black hole mass, stellar T 
and two dark halo parameters using the SAURON data and van 
der Marel data only, 3) we fit black hole mass, stellar T and 
two dark halo parameters using the globular cluster data only, 
and 4) we fit black hole mass, stellar T and two dark halo pa- 
rameters using the SAURON data, van der Marel data, and the 
globular cluster kinematics. For the last set of models we com- 
bine the constraints from the stellar and the GC kinematics by 
minimizing the sum Xstars + Xc;c- Note that because the stars 
and the GCs sample different spatial scales of M87, the best-fit 
models 2), 3) and 4) are not the same. 

For the models 1) without dark matter halo we have to deter- 
mine two free parameters: T andMeH (cf. equation[TJ. How- 
ever, when including a dark matter halo the amount of parame- 
ter coverage (with four parameters and two sets of potentials) is 
large. We start with a uniform, sparse grid in this four dimen- 
sional parameter space. Then we iteratively sample on a finer 
grid, focussing on the region(s) around the smallest % 2 . Specif- 
ically, during iteration n we determine all models closer to the 
best-fit than A\ 2 < 1. We then extend the sampling around 
these models by filling the 80 adjacent grid points around each 
of them. In this way, iteration n closes the surface in param- 
eter space around all models with A\ 2 < 1 of iteration n - 1 . 
After the last iteration our fine grid covers the entire A\ 2 < 1 
region. We also finely sampled a region that is not at the mini- 
mum to assure that the sampling is adequate in all areas of x 2 ; 
this denser sampling can be seen in the figures described below, 
and we find that our sampling procedure adequately represents 
the minimum \ 2 contours over the four parameters. 

We have run over 25 , 000 orbit libraries with about 25 , 000 
orbits per library, and each model takes approximately 1.5 
hours of cpu. We use the Texas Advanced Computing Center 
(TACC) at the University of Texas at Austin, where they have a 
cluster of over 5000 nodes. 

4. RESULTS 

We first discuss models for M87 excluding the contribu- 
tion from a dark halo. All black holes measured to date ex- 
clude the contribution from a dark halo (e.g., Gebhardt et al. 
2003), so these models provide a good comparison. The best- 
fitted MJL V is 10.2 ± 0.4, and the best fitted black hole mass is 
2.3(±0.6) x 10 9 Mq. Both of these values are consistent with 
previous dynamical determinations (as compiled in Kormendy 
& Richstone 1995, Cappellari et al. 2006). The model is com- 
pared to the data in Fig. 3 (dotted lines). It provides a good fit 
to the data. 

Including the dark halo significantly increases the parameter 
space that has to be explored. Given that we run over 10000 
models with four parameters it is difficult to visualize the % 2 
distribution with a single plot. Figure 4 plots the \ 2 versus 
each of the four parameters, including all values of the other 3 
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parameters. For each panel the contour that follows the enve- 
lope of minimum values is the marginalized \ 2 , and is what we 
use for estimation of the parameters and their uncertainties. 

4.1. x 2 Analysis 

The combined x 2 = Xstars + XGo which we use to constrain 
the gravitational potential. Figure 4 plots the \ 2 values for the 
stars only, globular clusters only and for the combined sample. 
The red line in the figure is for the stars (Xstars)' the blue line 
for the clusters Xgc)' the black line for the combination. The 
values on the y-axis are the sum of x 2 for the stars and clusters. 
In order to plot the individual \ 2 from the stars and clusters, 
we must vertically shift their values; since we only use an addi- 
tive shift, the relative \ 2 is preserved. A value of 26.3 is added 
to the x 2 f° r the stars, and 1 14.2 for that of the clusters. We 
note that the black line is not the sum of the blue and red lines, 
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since the model with the smallest x 2 in the sum will not be the 
same model for the one that gives the minimum when using 
either stars or clusters only. For the stars, with 52 LOSVDs 
that we use in the models and since each are generated from 
4 Gauss-Hermite parameters, the total number of parameters 
is 208. For the clusters, there are 11 LOSVDs, with about 4 
independent values, adding another 44 parameters. Thus, the 
reduced x 2 is about 0.65. This number is similar to what we 
find for orbit-based models (Gebhardt et al. 2003), which may 
be due to a correlation among the Gauss-Hermite coefficients 
(e.g., Houghton et al. 2006). 

It is clear from Figure 4 that the stars are primarily respon- 
sible for determining the black hole mass and that the clusters 
are responsible for determining the circular velocity of the dark 
matter. The stellar M/L comes primarily from the stars, but it 
has a significant degeneracy with the dark halo (see below). The 
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FIG. 4. — x 2 versus the four parameters. In each panel, we plot all models versus one of the parameters. A trace along the bottom ridge thus represents the 
marginalized \ 2 distribution, which we use to determine the best fit and uncertainties. The small points represent all models for the logarithmic potential. The solid 
lines are the marginalized x 2 values for three models. The red line is the \ 2 determined from stars only; the blue line comes from clusters only; the black line is 
from the combination. We note that the black line is not a sum of the red and blue lines. The value of the \ 2 (the vertical axes) are for the combined value, and that 
for stars and clusters have been shifted in order to see detail in the \ 2 contours. In the plot, for the stars, we add 26.3 to the x 2 , and for the clusters we add 1 14.2. 
The differential \ 2 values are preserved. 
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FIG. 5. — Plots of the four parameters against each other, designed to show possible degeneracies. The parameters are: M/L in the V-band, black hole mass in 
Mq , circular velocity of a logarithmic halo in km s , and core radius of the logarithmic halo in kpc. The plotting range is chosen to highlight the region around 
the best fit for each parameter. The small grey points represent the locations of all the models; there are many additional models beyond the plotting range that all 
provide significantly worse fits. The large black points are the locations of those models that are within the 68% confidence limit after marginalizing over the other 
parameters (i.e., A\ 2 = 1.0). The red points are those that are within the 95% confidence band. The only panel that shows an obvious correlation is the one with 
core radius against M/L (the upper-right plot), which is due to the trade between the mass of dark matter and stars. The correlation between stellar M/L and black 
hole mass is weak within the 68% bands shown here, but becomes significant for the 95% confidence band. 



scale radius of the dark halo is actually determined from both 
stars and clusters. 

Figure 4 provides the values of the best-fitted parameters 
and uncertainties. For the uncertainties we use A% 2 = 1 for 
each parameter when marginalizing over the other three param- 
eters. We determine the range of those models within the 68% 
limit and use the middle of the range as the best-fit and half 
of the range as the uncertainty. For black hole mass we find 
6.4(±0.5) x 10 9 M Q ; for stellar M/L v we get 6.3(±0.8); for 
logarithmic halo circular velocity we find 715(±15) km s , 
for dark halo scale radius we get 14.0(±2.0) kpc. Our best-fit 
model with dark matter halo is compared to the data in Figure 3 
(solid lines). In the outer parts, models with and without dark 
matter fit almost equally well to the SAURON data, reflecting 
that the information about the halo comes mostly from the GC 
kinematics. However, the model with dark matter halo fits bet- 



ter in the inner regions (both to the van der Marel and the inner 
SAURON data). This indicates that there is a significant change 
in the dynamical mass-to-light ratio fromM/Ly w 10 at r = 40" 
to M/Ly ~ 6 at r ~ 3". Towards the very center, the two models 
again become very similar, which is (1) due to the large uncer- 
tainties in the van der Marel data and (2) due to the degeneracy 
between stellar M/L and black-hole mass. 

Our results could be biased by the choice of the dark halo 
parameterisation. In particular, since we use a profile with a 
flat central density core, our models are entirely dominated by 
luminous mass in the inner regions. Then, in the center the 
luminous mass is degenerate with the black-hole mass and, in 
the outer parts, the luminous mass is degenerate with the halo 
(see also the discussion below), but there is no direct commu- 
nication between halo and black-hole. The situation would be 
different when using a cuspy halo profile, which could lead to 



8 



models with a significant amount of dark mass near the center, 
such that halo and black-hole mass could be directly degener- 
ate. To check for this, we have also run models for NFW. The 
grid used for the NFW models is rather coarse (we use about 
1/10 of the sampling used for the logarithmic potential) and a 
finer grid would be more informative, but the general shape of 
X 2 versus the four parameters is similar. It turns out that realis- 
tic NFW fits (that fit the GC kinematics) are not centrally con- 
centrated enough to dominate the mass in the inner regions and, 
thus, the mass profiles we get with NFW fits are very similar to 
those obtained with logarithmic halos. (Note also, that the cen- 
tral logarithmic slope of the deprojected light-profile is close to 
-1 and, thus, comparably steep as the cusp of the NFW core). 
That fits with NFW and logarithmic halos yield similar results 
has also been found in other early-type galaxies (Thomas et al. 
2005, 2007). Since the results are similar, we do not discuss 
results from NFW models. 

Figure 5 plots the relationships between the four parameters 
for the models. The parameter space is much larger than shown 
in this figure since we only focus on the region around the best 
fit. In each panel we include all models that are in that parame- 
ter range, irregardless of their \ 2 values. The color and symbol 
size, however, indicates those models that are within the 68% 
(Ax 2 = 1 after marginalizing over the other 3 parameters) and 
the 95% (Ax 2 = 4) confidence bands of the best-fitted model. 
The idea is to explore any correlations/degeneracies among the 
parameters. There is not an obvious correlation between M/L 
and black hole mass in the plotted region around the minimum 
X 2 . However, increasing the confidence band to 95% shows a 
correlation. This degeneracy is expected. Since the dynami- 
cal models only fit for the enclosed mass, as the contribution 
from the stars is increased in the central region, the contribu- 
tion from the black hole has to decrease. An obvious degener- 
acy within the 68% band is that between M/L and scale radius 
of the halo. As the scale radius is decreased, for a given cir- 
cular velocity, the contribution of the dark matter in the central 
regions is increased, thereby decreasing the contribution from 
the stars. These degeneracies are the essential reason why the 
dark halo properties affect the inferred black hole mass in M87. 
Beyond the plotting region, there are no points within the 68% 
confidence region, which allows us to set reasonable uncertain- 
ties on the parameters. 

4.2. Gravitational Potential 

Figure 6 plots the gravitational potential. The black lines 
are the potentials of the models that are within the 68% limit 
(Ax 2 = 1 for the marginalized parameters). We also plot the 
gravitational potential as inferred from the X-ray profile as pre- 
sented in Churazov et al. (2008). These have all been scaled to 
be zero at 330" as in Churazov et al. The red line is the model 
from Churazov et al. and it is clearly dissimilar from the best 
fit potential derived from the stellar and cluster dynamics. The 
green line is one of the models we used which is close to the 
gravitational potential from the X-rays; we use that model to 
represent the X-ray profile in subsequent analysis. Between the 
green line and the best-fitted model, Ax 2 = 18, which implies 
the potential derived from the X-ray is significantly poorer fit to 
the kinematics. We find a potential that is deeper than the X-ray 
derived potential. Churazov et al. (2008) explore a shock model 
for the X-rays in an attempt to explain the wiggles seen in the 
X-ray potential, but this model would not explain the large off- 
set. Additional galaxies with potentials derived from both stars 



and X-rays will be important to study. This difference in the 
potential may be specific to M87. 
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FIG. 6. — The gravitational potential versus radius for M87. The solid black 
lines are the dynamical models from this paper that are within the 68% confi- 
dence band of the best fit. The red noisy line is the gravitational potential as 
derived from the X-ray gas emission from Churazov et al. (2008). The green 
solid line is the parameterization from our models that well matches the X-ray 
profile (the parameters of the matched model are 6 X 10 9 Mq , 8.0, 800 km s , 
35 kpc for black hole mass, M/Ly, circular velocity and core radius, respec- 
tively. The difference in \ 2 between our representation of the X-ray potential 
and that of the best fit is 18. The units of the potential are as given in Churazov 
et al. and all models are scaled to R=330". 



4.3. Enclosed Mass and M/L Ratio 

Figure 7 plots the enclosed mass. The black lines represent 
the models that are within the 68% confidence limit. The en- 
closed mass flattens in the center due to the black hole and rises 
linearly due to the dark halo. The red line is the stellar contri- 
bution assuming M/L=5.0. The green line is the enclosed mass 
from our model which best matches the gravitational potential 
derived from the X-rays. In the potential derived from the X- 
rays (green line in Figure 6), there is a transition around 300" 
from the steep profile beyond that to a shallow profile inside of 
that radius. This transition manifests itself as a wiggle in the 
enclosed mass profile around 200-300". The enclosed mass 
derived from the stellar kinematics instead show a very smooth 
gradient from a few arcseconds out to 1000". If enclosed mass 
profiles are relatively smooth, as the one derived here, it will 
be fairly easy to constrain that profile using a non-parametric 
approach. This approach would be much cleaner to implement 
since the kinematics only measure enclosed mass. By trying to 
constrain the stellar contribution to the mass and the dark halo 
contribution leads to degeneracies as seen in the above Figures. 
However, this alternative approach would require some simpli- 
fication of the enclosed mass profile in order to use a reasonable 
amount of cpu time (e.g., the four parameters explored above 
required about 20,000 hours). 
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FIG. 7. — The mass profile for M87. The black lines represent the mod- 
els that are within the 68% confidence band of the best fit (as in Figure 6). 
The green line is the mass profile derived from our representation of the X-ray 
gravitational potential (i.e., the green line in Figure 6). The red line is the aver- 
age contribution from the stars, where we use the light profile in Figure 1 times 
6.3 (the best-fitted M/L). The mass profiles for the dynamical model show a 
smooth transition from 30 to 1000", whereas the X-ray profile shows a kink. 
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FIG. 8. — The integrated M/L versus radius. The red dotted lines represent 
the models that are within the 68% confidence band of the best-fitted potentials 
(black lines in Figure 6). The solid black lines are the same models excluding 
the contribution from the black hole. The green solid and dotted lines are our 
best representation of the X-ray potential without and with a black hole. The 
horizontal blue line is the M/L derived from stellar population models. From 
the red dotted lines, there is no radial region in M87 that has a constant M/L, 
since one always has to consider the contribution from either black hole or 
dark halo. 



Figure 8 plots the integrated mass-to-light ratio. We include 
all models within the 68% confidence limit of the best-fitted 
model. M/L profiles for these models are shown as the red dot- 
ted lines in the figure. The increase at small radius is from the 
black hole and the increase at large radius from the dark halo. 
There is no region in M87 where the M/L is constant. The solid 
black lines in Figure 8 are the M/L profiles including just the 
stars and dark halo, and are intended to show the range of stellar 
M/L values. Figure 8 also plots the M/L profile from our repre- 
sentation of the X-ray potential, both without (solid green) and 
with (dashed green) a black hole. The blue horizontal line is 
the best-fitted value of the stellar M/L from stellar population 
models (Cappellari et al. 2006). The stellar population models 
use Kroupa IMF with lower mass cut-off of 0.01. The stellar 
population models may have a large range of M/L values de- 
pending on the lower mass cut-off. The adopted value of 0.01 
is quite low and may bias the M/L's high if the lower cut-off is 
more massive. 

5. DISCUSSION 

There are a variety of results from this work, some deter- 
mined more robustly than others. Very robust results are 1) 
the black hole mass in M87 from stellar dynamical measures is 
over 2 times higher when including a dark halo, 2) the dynam- 
ical stellar M/L is 2x lower when including a dark halo, and 3) 
uncertainties on parameters significantly increase when includ- 
ing a dark halo. A model dependent result from this work, al- 
though still very significant, is that 4) the gravitational potential 
derived from X-rays for M87 is inconsistent with that derived 
from kinematics. Below we discuss each of these points and 
their implications. 

5.1. Black Hole Mass 

One of the more robust results from this work is the in- 
crease in the black hole due to inclusion of the dark matter 
halo. The reason is quite simple: dynamical models only mea- 
sure enclosed mass and if, in the central regions, one low- 
ers the stellar contribution with a smaller M/L, then that mass 
must be compensated by a corresponding increase in the black 
hole mass. For M87, since the light profile is so shallow, the 
dynamically-derived M/L value will contain information from 
stars that spend most of their time at large radii where the dark 
halo dominates. Thus, excluding a halo will result in pushing 
the M/L higher to match the observed velocity dispersions. One 
could partially compensate for this by excluding kinematics at 
large radii in a dynamical analysis. This may work somewhat 
for some galaxies, but as seen in Figure 8, there is no region in 
M87 where there is a constant M/L, causing a need to model 
simultaneously the black hole and the dark halo. 

A concern from this work could be that the dark halo param- 
eters are poorly constrained by the globular cluster kinematics, 
leading to a biased result due to the above degeneracy. We have 
therefore run a grid of models assuming the potential derived 
from the X-rays and fitting for black hole mass and stellar M/L. 
We find nearly identical results when fitting for the dark halo 
as well. The dark halo that we derived is only moderately more 
massive than the one from X-rays, but using either gives the 
same answer for the black hole mass. This is reassuring as it 
leads to a clean result for the black hole mass. Thus, we are con- 
fident that in M87 the black hole mass is 6.4(±0.5) x 10 9 M© 
derived from stellar kinematics. 

There are significant implications for having a black hole 
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mass this large. There has long been an issue as to why 
black hole masses derived from quasars can exceed 10 10 M@, 
whereas the masses derived from local galaxies never approach 
that. M87 and NGC4649 (Gebhardt et al. 2003, Humphrey 
et al. 2008) have the largest measured black holes around 
3 x 10 9 Mq. The space density of quasars with massive black 
holes imply that we should have many black holes approach- 
ing 10 10 out to about 100 Mpc (Shields et al. 2003). Many 
studies have pointed out this problem (e.g., Richstone et al. 
1998, Sheilds et al. 2003, Lauer et al. 2007, Bernardi et al. 
2007, Salviander et al. 2008). With a black hole mass of 
6.4 x 10 9 Mq, the agreement between the space density of the 
most massive black holes is in closer agreement with that de- 
rived from quasars. Because this change is due to an issue in 
the modeling (and not sample selection) and because only one 
object can have a such a large effect, it is not worthwhile at 
this point to do a detailed comparison. A much better approach 
would be to re-examine the other massive galaxies and study 
additional ones. 

Another implication of biased and larger black hole masses 
is the effect it will have on correlations of the black hole mass 
with global galaxy properties. The correlation of black hole 
with velocity dispersion (Gebhardt et al. 2000, Ferrarese & 
Merritt 2000) has led to a large amount of theoretical research 
as to the cause of the correlation and the value of the slope (e.g., 
Hopkins et al. 2007). If the measured black holes at the upper 
are biased low, this will obviously have a significant lever arm 
on the value of the slope and scatter. The black hole masses 
were already high and possibly hint at a curvature in the rela- 
tion (Wyithe 2006, and see the recent analysis of Gultekin et 
al. 2009). Correlations with other parameters would also need 
re-evaluating. Again, it is not prudent to re-analyse the current 
sample in regards to the correlations until the bias is understood 
well with both re-examinations and larger samples. 

There is a difference of the black hole mass measured with 
our analysis compared to the previous HST gas kinematic mea- 
sure. Macchetto et al. (1997) use gas kinematics measured with 
HST to derive a black hole mass of 3.2(±0.9) x 10 9 M . They 
use a distance of 15 Mpc. Transforming to our assumed dis- 
tance of 17.9 Mpc gives a mass of 3.8(±0.9) x 10 9 . The radii 
of the gas emission is in a region that is dominated by the black 
hole, so there is little difference if the stellar M/L changes, as 
we measure in this paper. For the black hole mass, 6.4 ± 0.5 
compared to 3.8 ±0.9 is only a 2-sigma difference and while 
that may not be statistically significant, it is important to con- 
sider possible causes for the difference. There have been very 
few studies of galaxies with black hole masses measured from 
both stars and gas. Of those that do have both measures, the 
results can vary greatly: IC1459 has the gas measure 8 times 
larger than that of the stars (Cappellari et al. 2002) and CenA 
has the stellar measure 7 times larger than that of the gas (Silge 
et al. 2005, Neumayer et al. 2007,), although Cappellari et 
al. (2009) show that they are consistent for CenA. The maser 
galaxy, NGC4258, has a consistent mass estimate (Siopis et al. 
2009), although using the masing gas is not the same as using 
the gas further out. Ho et al. (2002) discuss ways to optimize 
use of gas kinematics for measuring the black hole mass by fo- 
cussing on those galaxies with well-organized dust lanes. For 
M87, the gas disk is not very well defined (see Macchetto et al 
1997) and the uncertainties with the inclination can be large, but 
it would be difficult to have inclination explain all of the differ- 
ence. More significant concerns with using the gas kinematics 



is that the intrinsic density distribution is poorly constrained 
and the gas emission lines usually show large velocity broaden- 
ings that are not generally considered. For M87, Macchetto et 
al. include an estimate of the density distribution but there are 
clearly large uncertainties since many volume density distribu- 
tions can project to give similar projected densities. Macchetto 
et al. model the line centroids of the gas emission to determine 
the gravitational potential. For their assumed intrinsic gas dis- 
tribution (which has a hole in the center), they also match the 
large velocity dispersions seen in the emission lines; however, 
the unknown intrinsic distribution can strongly effect the gas 
centroids and dispersions. The problem is that it is not clear 
whether one should include the velocity broadening (thereby 
assuming a hot component to the gas orbits) or exclude them 
(assuming some type of turbulent motion), as discussed in Ver- 
does Kleijn et al. (2006). For M87, this effect could account 
for some of the difference. 

The resolution of the difference between the gas kinematics 
and the stellar kinematics will best come from stellar kinematic 
data taken with better spatial resolution. Whatever the value 
of the black hole mass however, the dynamical models show a 
strong effect between including and excluding a dark halo when 
the influence of the black hole is not well resolved. 

5.2. Stellar M/L 

The stellar M/Ly that we derive without a dark halo is 10 
(9.7 correcting for extinction) and with a dark halo is 6.3 (5.3 
correcting for extinction). This is a very significant difference, 
but it is obvious as to why that happens since the dynamical 
models simply are trading enclosed mass in stars with enclosed 
mass in dark matter. Adopting the gravitational potential as de- 
rived from the X-rays requires a stellar M/L of 8, but this model 
provides a significantly worse fit than our best fit. We compare 
to M/L as derived from stellar population models. Cappellari 
et al (2006) use Kroupa IMF to derive M/L 7 = 3.33 and find a 
dynamically-derived M/L/ = 6.1, using Ag = 0.1 and distance 
of 15.6 Mpc. Transforming to V-band and 17.9 distance im- 
plies M/Ly = 5.2 (stellar population) and M/Ly = 9.6 (dynami- 
cal). Our dynamically-derived M/L is in good agreement with 
that derived from the stellar population models. Also, when 
we exclude the dark halo, we find a very similar M/L as the 
dynamically-derived M/L from Cappellari et al., who do not 
include a dark halo. Thus, the M/L derived from stellar pop- 
ulation models supports our conclusion that the dark halo is 
important for deriving the proper M/L and, subsequently, the 
black hole mass. 

The physical understanding of the Fundamental Plane is af- 
fected by the dark matter fraction as a function of galaxy mass, 
which implies being able to measure the stellar M/L. Until more 
galaxies are examined at this level of detail (as in Thomas et al. 
2007b), we do not attempt to address implications for the Fun- 
damental Plane. 

5.3. Parameter Uncertainties 

As we increase in the number of parameters in the dynami- 
cal models, the uncertainties on those parameters will also in- 
crease. In fact, there are some parameters which we have not 
explored for M87. Two important ones include inclination and 
triaxiality. For inclination, if M87 is spherical, then it should 
not have a large effect. However, if M87 is very flattened 
and seen face-on, inclination will be important; although since 
galaxies with the mass of M87 are all seen as nearly round, 
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the likelihood is that they are spherical. In addition, M87 is 
flattened at large radii. Triaxiality is a bigger concern, since 
these large galaxies tend to be triaxial in numerical simulations 
(although no conclusive evidence yet either way from observa- 
tions). Triaxial models (van den Bosch et al. 2008) have only 
recently been developed, and it will be important to study M87 
as a triaxial system. The uncertainties play a key role in un- 
derstanding the physical nature of the black hole correlations. 
Including a dark halo already increased the uncertainties on the 
black hole and particularly on the stellar M/L. Thus, until ad- 
ditional analysis on the degeneracies is explored and a larger 
sample is included, it would be prudent to include a modest 
systematic bias uncertainty for the black hole mass estimates 
when studying the correlations and their intrinsic scatter. 

5.4. Dark Halo and X-ray Properties 

It is generally assumed that the X-ray emitting gas can be 
used as a robust tracer of the gravitational potential, and there 
have been studies that show rough agreement between the two 
(as shown in Churazov et al. 2008). However, all of the dy- 
namical studies rely on fairly restrictive assumptions about the 
stellar distribution function: Wu & Tremaine (2006) assume a 
two-integral distribution function, McLaughlin (1999) assume 
spherical isotropic distribution function, and both Kronawitter 
et al. (2000) and Romanowsky & Kochanek (2001) assume 
a spherical distribution. Given that the stellar density profile 
at large radii is very flattened and that the globular clusters 
show significant rotation, it is important to explore more gen- 
eral models as the ones presented here. We find rough agree- 
ment with the enclosed mass at some locations in M87, but the 
shape of the potential profile is quite different. The gravita- 
tional potential in Figure 6, in fact, shows fairly significant dif- 
ferences between the two derived potentials. This figure, how- 
ever, is not as fair a comparison since the dynamical models 
directly measure the enclosed mass and it is more instructive 
to focus on that aspect. Since the enclosed mass is a derivative 
of the potential, one needs to focus on the change in shape of 
Figure 6; for R < 40", the derivative of the potential is very 
similar, which subsequently leads to good agreement with the 
enclosed mass. At R > 600", there is also good agreement in 
the enclosed mass (and derivative of the potential). But it is the 
radii from 40 < R < 600 where there is the strongest disagree- 
ment, since the potential derived from X-rays shows a strong 
kink there while that derived from stellar kinematics shows a 
smooth transition. This region is where there is no stellar data 
with modest coverage by globular cluster kinematics; however, 
given that we are using parametric models for the dark matter 
profile, it is a robust result that the X-ray potential is a sig- 
nificantly worse fit to the data than the best-fitted logarithmic 
potential. It is clear though that additional data at these radii 
are very much needed. We note that the kinematics of Sembach 
& Tonry (1996) are more consistent with our best-fitted profile 
than with the X-ray profile; however, the uncertainties are large 
and better data is needed to make the comparison conclusive. 

The dark halo parameters are, for a logarithmic potential, 
a circular velocity of 715 ± 15 km s"'and a core radius of 
14 ± 2 kpc. Kronawitter et al. (2000) finr a similar halo for 
M87. This halo is also similar to those measured in Thomas 
et al. (2007b) and Gerhard et al. (2001) in terms of the cir- 
cular velocity but M87 has a smaller core radius for galaxies 
of somewhat smaller mass. We also ran a smaller grid of po- 
tentials that follow an NFW profile. The best fitted parameters 



were around a concentration of 15 and a scale radius of 50 kpc. 
We do not provide uncertainties since we did not explore a very 
dense grid. Using either NFW or logarithmic, the potential de- 
rived from X-rays was statistically excluded. 

M87 may not be the best case where X-ray gas emission 
provides an adequate representation of the potential. The X- 
ray potential shows significant deviations from a smooth curve, 
and Churazov et al. (2008) use these deviations to derive a jet 
model, with cooling and heating regions. This model, however, 
does not explain the difference at 40 < R < 600". It would be 
important to study other galaxies that do not show the strong 
deviations seen in M87, but it is possible that the X-ray gas in 
M87 may not be an adequate tracer of the potential. 

5.5. Next Steps 

There are a variety of ways to move forward. On the data 
side, there is a strong need for observations at small radii using 
observations assisted with adaptive optics on a large telescope. 
The faint central surface brightness of M87 made it not possible 
to observe with the spectrographs on HST, but it is an ideal tar- 
get for ground-based telescopes. It should be straight-forward 
to detect a black hole with as high a mass as presented in this 
paper. At larger radii, there is a critical need to obtain stellar 
kinematics in the region from 40 < R < 300". This is feasible 
using current ground-based instruments (e.g., VIRUS-P on the 
2.7 m telescope at McDonald Observatory, Hill et al. 2008). 
The large radial data will strongly limit dark halo models, and 
allow us to not rely as heavily on the kinematics from the small 
number of globular clusters. Both of these observational av- 
enues, at small and large radii, are in progress. 

On the modeling side, it would be worthwhile to include 
more general models, such as triaxiality, and an exploration of 
inclination. For this paper, we use a parametric form for the 
dark halo and there is little theoretical justification for such an 
approach. A cleaner approach would be to use a non-parametric 
estimate for the enclosed mass, such as a mass profile defined 
in radial bins. This method would directly include any change 
of mass-to-light ratio of the stars. The enclosed mass is directly 
measured from the kinematics, and then interpretation of how 
to split between stars and dark halo could come subsequently. 
Given that the enclosed mass profile is relatively smooth (Fig- 
ure 7), it should be an easy matter to constrain it. However, 
with just 4 parameters, it was already a large computer effort, 
which could obviously increase dramatically as the number of 
fitted mass bins increase. 

There is a strong need to explore additional galaxies. It 
could be that since M87 has no region where there is a con- 
stant M/L (see Figure 8), that the degeneracies between black 
hole and dark halo, mediated by the stellar M/L become very 
large. Galaxies with smaller black holes and/or smaller rela- 
tive fraction of dark halo should be studied with this type of 
analysis. 

Whether or not the dark halo parameters presented here cor- 
rectly reflect reality (we obviously agree that they do), we per- 
formed a simple test of measuring the black hole with the same 
modeling code including and excluding a dark halo. In a differ- 
ential sense, the importance of including the dark halo is one of 
the strongest and cleanest results. All black holes measured, to 
date, including those of the authors (Gebhardt et al. 2003) ex- 
clude a dark halo (but see the discussion in Nowak et al. 2008). 
While we expect that M87 shows the largest effect due to its 
shallow light profile, we are at the point now where the details 
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of the black hole correlations are important and understanding 
all biases are therefore essential. 
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